; docformat = 'rst'
;
;+
;
; :Purpose:
;   Simple Bayesin inversion based on given measurements, uncertainties and sensitivites
;   
; :Inputs:
;   sim_name: (string) simulation name
;   
; :Requires:
;   Combined sensitivity from *cels_sensitivity.pro*, stored in {input_directory}/{sim_name}/sensitivity.sav
;   
;   Measurements in {input_directory}/{sim_name}/measurements/measurements.sav
;   
;   State definitions in {input_directory}/{sim_name}/state.sav
;   
;   Region definitions in {input_directory}/{sim_name}/regions.sav
;   
; :Outputs:
;   {input_directory}/{sim_name}/optimized.sav contains optimized state vector (x) and optimized observations (y_opt) 
;   
; :History:
; 	Written by: Matt Rigby, MIT, Aug 19, 2011
;
;-
pro cels_invert, sim_name

  compile_opt idl2
  on_error, 2
  
  print, "CELS_INVERT: Bayesian inversion"

  StartY=cels_get_parameter(sim_name, 'STARTY')
  EndY=cels_get_parameter(sim_name, 'ENDY')
  IC_uncertainty=cels_get_parameter(sim_name, 'IC_uncertainty')
  NLR_uncertainty=cels_get_parameter(sim_name, 'NLR_uncertainty')
  LR_uncertainty=cels_get_parameter(sim_name, 'LR_uncertainty')

  restore, cels_filestr(/Input, sim_name + '/state.sav')
  restore, cels_filestr(/Input, sim_name + '/regions.sav')
  restore, cels_filestr(/Input, sim_name + '/sensitivity.sav')

  ;DEFINE UNCERTAINTIES
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    
  P_apd=(lr_uncertainty*x_ap)^2
  p_apd[0:n_IC-1]=(IC_uncertainty*x_ap[0:n_IC-1])^2
  p_apd[n_IC:N_NLRT-1]=(NLR_uncertainty*x_ap[n_IC:n_NLRT-1])^2

  p1_apd=1./p_apd
  p1=diag_matrix(P1_apd)
  
  R1=diag_matrix(1./(y_error^2))

  ;INVERT
  P=invert(transpose(H)##R1##H + P1)
  x=x_ap + P##transpose(H)##R1##(y-H##x_ap)

  caldat, x_time, dummy, dummy, x_year
  
  for yi=StartY, EndY-1 do begin
    print, 'CELS_INVERT: ', yi, ' emissions: ', total(x[where(x_year eq yi)]*365.*24.*3600.)/1.e6, ' Gg/yr', $
      format='(A, I4, A, f6.2, A)'
  endfor

  y_opt=h##x

  save, filename=cels_filestr(/Input, sim_name + '/optimized.sav'), $
    x, y_opt

  
end